#############################
# Carsten Minh Vo
# Geopolitical Competition's effect on State Engagement with the Global Political Economy
# 2023
############################

#################################
# Reading data and data management
#################################

rm(list = ls())


library(foreign)
library(plm)
library(ggplot2)
library(dplyr)
library(stargazer)
library(stringr)
library(tidyr)
library(MASS)
library(gridExtra)
library(tidyverse)
library(readr)
library(readxl)
library(lmtest)
library(sandwich)
library(car)
library(pscl)
library(AER)
library(lfe)



###################################
# Uploading Data and filtering data
###################################

MyData <- read.csv("MyData.csv", stringsAsFactors = F)
names(MyData)



# Converting my data
# Create a vector with the desired factor labels
factor_labels <- c("Autocracy", "Closed Anocracy", "Open Anocracy", "Democracy", "Full Democracy")
# Convert the "polity2" variable into a factor with the specified labels
MyData$polity2 <- cut(MyData$polity2,
                      breaks = c(-Inf, -5, 0, 5, 9, Inf),
                      labels = factor_labels,
                      include.lowest = TRUE)

# Check the updated factor levels
levels(MyData$polity2)
# Check the updated factor levels
levels(MyData$polity2)
MyData$polity2 <- as.factor(MyData$polity2)
MyData$Income <- as.factor(MyData$Income)


# Filter my data
MyData <- subset(MyData, Military_Spending < 50)
MyDataFDIIn <- subset(MyData, FDI_Inflow > 0.00000)
MyDataFDIOut <- subset(MyData, FDI_Outflow > 0.00000)




# Create new variables with the natural logarithm
MyData$Log_Geopolitical_Competition <- log(MyData$Geopolitical_Competition)
MyDataFDIIn$Log_Geopolitical_Competition <- log(MyDataFDIIn$Geopolitical_Competition)
MyDataFDIOut$Log_Geopolitical_Competition <- log(MyDataFDIOut$Geopolitical_Competition)

MyData$Log_Military_Spending <- log(MyData$Military_Spending)
MyDataFDIIn$Log_Military_Spending <- log(MyDataFDIIn$Military_Spending)
MyDataFDIOut$Log_Military_Spending <- log(MyDataFDIOut$Military_Spending)

MyData$Log_Trade <- log(MyData$Trade)

MyDataFDIIn$Log_FDI_Inflow <- log(MyDataFDIIn$FDI_Inflow)
MyDataFDIOut$Log_FDI_Outflow <- log(MyDataFDIOut$FDI_Outflow)

MyData$Log_Population <- log(MyData$Population)
MyDataFDIIn$Log_Population <- log(MyDataFDIIn$Population)
MyDataFDIOut$Log_Population <- log(MyDataFDIOut$Population)

MyData$Log_GDP <- log(MyData$GDP)
MyDataFDIIn$Log_GDP <- log(MyDataFDIIn$GDP)
MyDataFDIOut$Log_GDP <- log(MyDataFDIOut$GDP)

MyData$Log_Rivals <- log(MyData$Rivals)
MyDataFDIIn$Log_Rivals <- log(MyDataFDIIn$Rivals)
MyDataFDIOut$Log_Rivals <- log(MyDataFDIOut$Rivals)




######################################################################
########## Regression Analysis
######################################################################
### Military Spending
###################################
# Fit multiple regression models
model_Log_Military_trade <- lm(Log_Trade ~ Log_Military_Spending + polity2 + Income + Log_Population, data = MyData)
model_Log_Military_fdi_inflow <- lm(Log_FDI_Inflow ~ Log_Military_Spending + polity2 + Income + Log_Population, data = MyDataFDIIn)
model_Log_Military_fdi_outflow <- lm(Log_FDI_Outflow ~ Log_Military_Spending + polity2 + Income + Log_Population, data = MyDataFDIOut)


bptest(model_Log_Military_trade)
coeftest(model_Log_Military_trade, cvov = vcovHC)
bptest(model_Log_Military_fdi_inflow)
coeftest(model_Log_Military_fdi_inflow, cvov = vcovHC)
bptest(model_Log_Military_fdi_outflow)
coeftest(model_Log_Military_fdi_outflow, cvov = vcovHC)


# View regression results
summary(model_Log_Military_trade)
summary(model_Log_Military_fdi_inflow)
summary(model_Log_Military_fdi_outflow)



###################################
# Dianostic testing
###################################
# Residual analysis for the "Military_Trade" regression model
par(mfrow = c(2, 2)) # Arrange plots in a 2x2 grid

plot(model_Log_Military_trade, which = 1)  # Residuals vs. Fitted Values
plot(model_Log_Military_trade, which = 2)  # Normal Q-Q plot
plot(model_Log_Military_trade, which = 3)  # Scale-Location plot
plot(model_Log_Military_trade, which = 5)  # Residuals vs. Leverage

# Residual analysis for the "Military_FDI_inflow" regression model

plot(model_Log_Military_fdi_inflow, which = 1)  # Residuals vs. Fitted Values
plot(model_Log_Military_fdi_inflow, which = 2)  # Normal Q-Q plot
plot(model_Log_Military_fdi_inflow, which = 3)  # Scale-Location plot
plot(model_Log_Military_fdi_inflow, which = 5)  # Residuals vs. Leverage

# Residual analysis for the "Military_FDI_Outflow" regression model

plot(model_Log_Military_fdi_outflow, which = 1)  # Residuals vs. Fitted Values
plot(model_Log_Military_fdi_outflow, which = 2)  # Normal Q-Q plot
plot(model_Log_Military_fdi_outflow, which = 3)  # Scale-Location plot
plot(model_Log_Military_fdi_outflow, which = 5)  # Residuals vs. Leverage



# Calculate robust standard errors clustered by country
robust_se_model_Log_Military_trade <- sqrt(diag(vcovCL(model_Log_Military_trade, cluster = MyData$country)))
robust_se_model_Log_Military_fdi_inflow <- sqrt(diag(vcovCL(model_Log_Military_fdi_inflow, cluster = MyDataFDIIn$country)))
robust_se_model_Log_Military_fdi_outflow <- sqrt(diag(vcovCL(model_Log_Military_fdi_outflow, cluster = MyDataFDIOut$country)))

# Output the regression results with robust standard errors
################################################################
####[RESULTS FOR MODEL 1.1, 2.1, AND 3.1, SECTION 6.4]
################################################################
summary(model_Log_Military_trade, robust = TRUE, robust.se = robust_se_model_Log_Military_trade)
summary(model_Log_Military_fdi_inflow, robust = TRUE, robust.se = robust_se_model_Log_Military_fdi_inflow)
summary(model_Log_Military_fdi_outflow, robust = TRUE, robust.se = robust_se_model_Log_Military_fdi_outflow)
################################################################

# Calculate VIF for the regression model
vif_values_Military_Trade <- vif(model_Log_Military_trade)
vif_values_Military_FDIIn <- vif(model_Log_Military_fdi_inflow)
vif_values_Military_FDIOut <- vif(model_Log_Military_fdi_outflow)


# View the VIF values
print(vif_values_Military_Trade)
print(vif_values_Military_FDIIn)
print(vif_values_Military_FDIOut)


# Number of obs
N1.1 <- nobs(model_Log_Military_trade)
N1.2 <- nobs(model_Log_Military_fdi_inflow)
N1.3 <- nobs(model_Log_Military_fdi_outflow)

print(N1.1)
print(N1.2)
print(N1.3)



###################################
###### Regression Analysis - Geopolitical Competition
###################################
###################################
# Log Data
###################################
# Fit multiple regression models
model_LogGC_trade <- lm(Log_Trade ~ Log_Geopolitical_Competition + Income + Log_Population, data = MyData)
model_LogGC_fdi_inflow <- lm(Log_FDI_Inflow ~ Log_Geopolitical_Competition + Income + Log_Population, data = MyDataFDIIn)
model_LogGC_fdi_outflow <- lm(Log_FDI_Outflow ~ Log_Geopolitical_Competition + Income + Log_Population, data = MyDataFDIOut)




# View regression results
summary(model_LogGC_trade)
summary(model_LogGC_fdi_inflow)
summary(model_LogGC_fdi_outflow)


###################################
# Dianostic test
###################################
# Residual analysis for the "GC_Trade" regression model
par(mfrow = c(2, 2)) # Arrange plots in a 2x2 grid

plot(model_LogGC_trade, which = 1)  # Residuals vs. Fitted Values
plot(model_LogGC_trade, which = 2)  # Normal Q-Q plot
plot(model_LogGC_trade, which = 3)  # Scale-Location plot
plot(model_LogGC_trade, which = 5)  # Residuals vs. Leverage

# Residual analysis for the "GC_FDI_inflow" regression model

plot(model_LogGC_fdi_inflow, which = 1)  # Residuals vs. Fitted Values
plot(model_LogGC_fdi_inflow, which = 2)  # Normal Q-Q plot
plot(model_LogGC_fdi_inflow, which = 3)  # Scale-Location plot
plot(model_LogGC_fdi_inflow, which = 5)  # Residuals vs. Leverage

# Residual analysis for the "GC_FDI_Outflow" regression model

plot(model_LogGC_fdi_outflow, which = 1)  # Residuals vs. Fitted Values
plot(model_LogGC_fdi_outflow, which = 2)  # Normal Q-Q plot
plot(model_LogGC_fdi_outflow, which = 3)  # Scale-Location plot
plot(model_LogGC_fdi_outflow, which = 5)  # Residuals vs. Leverage



# Calculate robust standard errors clustered by country
robust_se_model_LogGC_trade <- sqrt(diag(vcovCL(model_LogGC_trade, cluster = MyData$country)))
robust_se_model_LogGC_fdi_inflow <- sqrt(diag(vcovCL(model_LogGC_fdi_inflow, cluster = MyDataFDIIn$country)))
robust_se_model_model_LogGC_fdi_outflow <- sqrt(diag(vcovCL(model_LogGC_fdi_outflow, cluster = MyDataFDIOut$country)))

################################################################
# Output the regression results with robust standard errors
#### [RESULTS MODEL 1.3, 2.3, AND 3.3, APPENDIX A.2]
################################################################
summary(model_LogGC_trade, robust = TRUE, robust.se = robust_se_model_LogGC_trade)
summary(model_LogGC_fdi_inflow, robust = TRUE, robust.se = robust_se_model_LogGC_fdi_inflow)
summary(model_LogGC_fdi_outflow, robust = TRUE, robust.se = robust_se_model_model_LogGC_fdi_outflow)


# Calculate VIF for the regression model
vif_values_GC_Trade <- vif(model_LogGC_trade)
vif_values_GC_FDIIn <- vif(model_LogGC_fdi_inflow)
vif_values_GC_FDIOut <- vif(model_LogGC_fdi_outflow)


# View the VIF values
print(vif_values_GC_Trade)
print(vif_values_GC_FDIIn)
print(vif_values_GC_FDIOut)



# Number of obs
N1.1 <- nobs(model_LogGC_trade)
N1.2 <- nobs(model_LogGC_fdi_inflow)
N1.3 <- nobs(model_LogGC_fdi_outflow)

print(N1.1)
print(N1.2)
print(N1.3)


###################################
# Rivals - Multiple Linear Regression
###################################

model_Rivals_trade <- lm(Log_Trade ~ Rivals + polity2 + Income + Log_Population, data = MyData)
model_Rivals_fdi_inflow <- lm(Log_FDI_Inflow ~ Rivals + polity2 + Income + Log_Population, data = MyDataFDIIn)
model_Rivals_fdi_outflow <- lm(Log_FDI_Outflow ~ Rivals + polity2 + Income + Log_Population, data = MyDataFDIOut)


# View regression results
summary(model_Rivals_trade)
summary(model_Rivals_fdi_inflow)
summary(model_Rivals_fdi_outflow)

# Diagnostics

plot(model_Rivals_trade, which = 1)  # Residuals vs. Fitted Values
plot(model_Rivals_trade, which = 2)  # Normal Q-Q plot
plot(model_Rivals_trade, which = 3)  # Scale-Location plot
plot(model_Rivals_trade, which = 5)  # Residuals vs. Leverage

plot(model_Rivals_fdi_inflow, which = 1)  # Residuals vs. Fitted Values
plot(model_Rivals_fdi_inflow, which = 2)  # Normal Q-Q plot
plot(model_Rivals_fdi_inflow, which = 3)  # Scale-Location plot
plot(model_Rivals_fdi_inflow, which = 5)  # Residuals vs. Leverage

plot(model_Rivals_fdi_outflow, which = 1)  # Residuals vs. Fitted Values
plot(model_Rivals_fdi_outflow, which = 2)  # Normal Q-Q plot
plot(model_Rivals_fdi_outflow, which = 3)  # Scale-Location plot
plot(model_Rivals_fdi_outflow, which = 5)  # Residuals vs. Leverage



# BP-test
# Perform the Breusch-Pagan test
bp_test_model_Rivals_trade <- bptest(model_Rivals_trade)
bp_test_model_Rivals_fdi_inflow <- bptest(model_Rivals_fdi_inflow)
bp_test_model_Rivals_fdi_outflow <- bptest(model_Rivals_fdi_outflow)


# Print the test results
print(bp_test_model_Rivals_fdi_outflow)
print(bp_test_model_Rivals_fdi_inflow)
print(bp_test_model_Rivals_trade)

# Calculate robust standard errors clustered by country
robust_se_model_Rivals_trade <- sqrt(diag(vcovCL(model_Rivals_trade, cluster = MyData$country)))
robust_se_model_Rivals_fdi_inflow <- sqrt(diag(vcovCL(model_Rivals_fdi_inflow, cluster = MyDataFDIIn$country)))
robust_se_model_Rivals_fdi_outflow <- sqrt(diag(vcovCL(model_Rivals_fdi_outflow, cluster = MyDataFDIOut$country)))

################################################################
# Output the regression results with robust standard errors
#### [RESULTS FOR MODEL 1.2, 2.2, AND 3.2, SECTION 6.4]
################################################################
summary(model_Rivals_trade, robust = TRUE, robust.se = robust_se_model_Rivals_trade)
summary(model_Rivals_fdi_inflow, robust = TRUE, robust.se = robust_se_model_Rivals_fdi_inflow)
summary(model_Rivals_fdi_outflow, robust = TRUE, robust.se = robust_se_model_Rivals_fdi_outflow)


model_Rivals_trade <- lm(Log_Trade ~ Rivals + polity2 + Income + Log_Population, data = MyData)
robust_se_model_Rivals_trade <- sqrt(diag(vcovCL(model_Rivals_trade, cluster = MyData$country)))
summary(model_Rivals_trade, robust = TRUE, robust.se = robust_se_model_Rivals_trade)


# Calculate VIF for the regression model
vif_values_Rivals_Trade <- vif(model_Rivals_trade)
vif_values_Rivals_FDIIn <- vif(model_Rivals_fdi_inflow)
vif_values_Rivals_FDIOut <- vif(model_Rivals_fdi_outflow)


# View the VIF values
print(vif_values_Rivals_Trade)
print(vif_values_Rivals_FDIIn)
print(vif_values_Rivals_FDIOut)

# Number of obs
N2.1 <- nobs(model_Rivals_trade)
N2.2 <- nobs(model_Rivals_fdi_inflow)
N2.3 <- nobs(model_Rivals_fdi_outflow)

print(N2.1)
print(N2.2)
print(N2.3)


###################################



